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Abstract 



The time dependence of the dark energy density can be an important clue 
to the nature of dark energy in the universe. We show that future supernova 
data from dedicated telescopes (such as SNAP), when combined with data 
of nearby supernovae, can be used to determine how the dark energy density 
px{z) depends on redshift, if px{z) is not too close to a constant. For quan- 
titative comparison, we have done an extensive study of a number of dark 
energy models. Based on these models we have simulated data sets in order 
to show that we can indeed reconstruct the correct sign of the time depen- 
dence of the dark energy density, outside of a degeneracy region centered on 
1 + wq = —wiZmax/'i (where z^ax is the maximum redshift of the survey, 
e.g., Zmax = 1-7 for SNAP). We emphasize that, given the same data, one can 
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obtain much more information about the dark energy density directly (and 
its time dependence) than about its equation of state. 

PACS numbers: 98.80.-k; 98.80.Es; 97.60.Bw 



Typeset using REVT^ 



I. INTRODUCTION 



Most of the energy in our universe is of unknown nature to us. The amount of this dark 
energy has been determined by recent experiments, including the Wilkinson Microwave 
Anisotropy Probe (WMAP) sateUite observations [1] of the anisotropy in the cosmic mi- 
crowave background radiation. Our universe is spatially flat (the three-dimensional equiva- 
lent of a two-dimensional plane), with roughly 27% matter and 73% dark energy. Determin- 
ing the nature of this dark energy is one of the major fundamental challenges in astronomy 
and physics today. 

There are many plausible candidates for dark energy. For example, (1) a cosmolog- 
ical constant, i.e., constant vacTium energy originally proposed by Einstein in his equa- 
tions of general relativity, (2) a time dependent vacuum energy, or scalar field known as 
"quintessence", that evolves dynamically with time [2], or (3) modified Friedmann equa- 
tion, e.g. the Cardassian models [3-5], that could result as a consequence of our observable 
universe living as a 3-dimensional brane in a higher dimensional universe. Other proposed 
modifications to the Friedmann equation include [6]. The time-dependence of the density of 
dark energy can reveal the nature of dark energy at a fundamental level. 

A powerful probe of dark energy is type la supernovae (SNe la), which can be used as 
cosmological standard candles to measure how distance depends on redshift in our universe. 
Observations of SNe la have revealed the existence of dark energy in the universe [7,8]. 
Current SN la data are not yet very constraining on the nature of the dark energy [24] . 

The distance-redshift relation of observed supernovae depends on the nature of dark 
energy. Most researchers have chosen to parametrize dark energy by its equation of state 
parameter. However, it has been shown [9,10] that it is extremely difficult to constrain 
the time-dependence of the dark energy equation of state using supernova searches (or any 
other technique relying on the luminosity distance); hence one might worry that one cannot 
differentiate between different dark energy models. Fortunately, it has been shown that one 
can do much better if one parametrizes the dark energy by its density directly, instead of its 
equation of state [11-14]. The dark energy density is a more fundamental parameter than the 
dark energy equation of state parameter. Obtaining the equation of state parameter requires 
one to perform an additional integral (compared to obtaining the dark energy density); this 
integral smears out much of the information one could otherwise learn. Hence, given the 
same data, the uncertainties of the constraints on the dark energy density should be smaller 
than that of the constraints on the dark energy equation of state. 

In this paper we focus on extracting information about the dark energy density directly. 
We will show how well, using future supernova data, one can determine whether the dark 
energy density changes with time, and whether it increases or decreases with time. 

We begin in Section II with the basic equations for using supernovae to study dark energy. 
In Section III we present four theoretical models which we will study to see how well we can 
reconstruct the time dependence of the dark energy: Models 1, 2, and 3 have dark energy 
density that is constant, increasing, and decreasing in time respectively. As our fourth set 
of models we consider those parametrized by an equation of state Wx{z) = Wq + WiZ. In 
Section IV, we simulated SNIa data for these models. In Section V, we use the adaptive 
iteration method to see how well we can reconstruct the time dependence of the dark energy 
density for these models: we use three test functions with different time dependences to see 
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which one best matches the data. For each model we then run 1000 Monte Carlo samples 
to obtain error bars for our fit. The results are presented in Section VI, followed by the 
conclusions. 



II. BASIC EQUATIONS 

Type la supernovae (SNe la) are our best candidates for cosmological standard candles, 
because they can be calibrated to have small scatters in their peak luminosity [15,16]. The 
distance modulus for a standard candle at redshift z is 

ix^{z) = m - M = 5 log j + 25, (1) 

where m and M are the apparent and absolute magnitudes of the standard candle, and 
di,{z) is its luminosity distance. 

In a flat Friedmann-Robertson- Walker universe (which we assume in this paper since it 
is strongly suggested by current CMB data [17,1]) the luminosity distance diXz) is given by 
[18] 

(l + z)c f dz' 

^^^"^ = Jo wr ^ ^ 

where Hq is the present value of the Hubble constant, and 

E{z) = ^n„(l + ^)3 + 0^/^(^). (3) 

Here Jl^ is the present value of the matter density in unit of the critical density pcrit = 
?>Hl/{%T^G), and Qx is the present value of the dark energy density in the same units. 
The condition for a flat universe imposes the relation 

Qm + ^X^ 1. (4) 

The dark energy density function 

h(^) = ^ (5) 

describes the redshift dependence of the dark energy density px{z). 

We note that Cardassian models [3-5] contain matter and radiation only (no vacuum), 
so that in those models, Jl^ and Q^x are used to refer to effective observed matter density 
and dark energy density respectively. 

For a given cosmological model with dark energy density px{z), or dark energy function 
fx{z), we can compare the measured distance modulus of SNe la at various redshifts with the 
predicted distance modulus of a standard candle at these redshifts. A systematic comparison 
spanning all plausible models yields constraints on the dark energy density px{z)- 
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III. DARK ENERGY DENSITY FUNCTIONS 



In this paper, we are interested in finding what information future SN la data can 
give us about the redshift dependence of the dark energy density. For this purpose, we 
consider four classes of dark energy densities: (1) constant with redshift, (2) increasing with 
redshift, (3) decreasing with redshift, and (4) a grid of models which includes some that 
are non-monotonic with respect to redshift. For each of these four classes we choose simple 
representative models. Some of the models we choose can be parametrized by a simple 
equation of state. 



Wx{z) ^ Wo + WiZ 



(6) 



in order to allow simple comparison with previous work in the literature. The simple 

paramctrization of Eq.(6) also allows us to estimate how small a deviation from a con- 
stant dark energy density can be determined by our technique. The dark energy density 
is constant only for wq = —1 and Wi = 0; any other values of wq or of wi parametrize 
deviations from a constant. However, we stress that it is the dark energy density itself, 
Px{z), that we extract from simulated data, as it is the more fundamental and more easily 
extracted quantity. 



Four Models 

The four sample theoretical models we consider are (see Table 1): 

Table 1 
Dark Energy Models 



Model 


model parameters 




Model 1: ACDM 
Model 2: quintessence model 
Model 3: MP Cardassian model 
Model 4: grid of models 


rirn ^ 0.3, flA = 0.7 

= 0.3, Wq(z) = — 1 + 0.5 Z 

0^* = 0.3, n = 0.2, q = 2 

wx{z) = Wo + WiZ with arbitrary wq, Wi 


p'x{z) = 
P'xi^) < 



(1) MODEL 1: For a constant dark energy density, p'x{z) = 0, we have a cosmological 
constant model {wq — —1 and Wi = 0). 

(2) MODEL 2: For an increasing dark energy density, p'x{z) > 0, we choose a 
quintessence model with equation of state wx{z) = — 1 + 0.5^;. Popular quintessence models 
in which the quintessence field tracks the matter field have p'x{z) > [10]. 

(3) MODEL 3: For a decreasing dark energy density, p'x{z) < 0, we choose a Modified 
Polytropic (MP) Cardassian model [4] with n — 0.2 and q = 2, 



Px{z) = Pcrit^m(l + 



1 + 



'I + zfq{l-n) 



- 1 



(7) 



MP Cardassian models can have either p'x{z) > or p'x{z) < 0. Our previous paper [5] 
shows the regions of parameter space that fall into the two regimes. We also discussed there 
that MP Cardassian models can be found with wx — Px/ px < —1 but with w — p/ p > — 1, 
so that the dominant energy condition holds (here, w refers to the total energy density 
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whereas Wx refers only to the new component in the Friedmann equation that mimics a 
dark energy). An effective wx < — 1 is consistent with recent CMB and large scale structure 
data [19,20]. 

(4) MODEL 4: We consider a grid of models of the form Eq.(6), for -1.2 < Wq < -0.5, 
and —1.5 < wi < 0.5, and a grid spacing of Awq — 0.1 (0.2 for wq < —1) and Awi — 
0.1. This represents a total of 147 models. This grid includes models with p'x{z) that is 
monotonically increasing or decreasing with redshift, as well as models with dark energy 
density that is non-monotonic with redshift. In a moment we will show which values of Wq 
and Wi correspond to a non-monotonic dark energy density. 

Although models of the form of Eq.(6) do not correspond exactly to physically motivated 
models, it is interesting to note that they can approximate a wide range of models. For 
example, Model 3 (MP Cardassian model with n = 0.2, q = 2. and Qm = 0.3) can be 
roughly approximated by a dark energy model with ■u;x(^) = —1.10—0.352;. The approximate 
equivalence of these models does not extend to the behavior of dark energy fluctuations but 
is limited to the average properties of the energy density. 

The dark energy density corresponding to Eq.(6) is [11] 

Px{z) = px{0)e^'"'' (1 + 2)3(^+"'°-^i). (8) 
Its derivative with respect to redshift follows as 

/ / X / 3{l+Wo + Wiz) 

Pxi^)^ Px [z) — . (9) 

i + z 

Thus p'x{z) has the same sign as 1 + + Wiz. 

Nonmononotic models: Since z must be positive, models with (i) wq < —1 and wi > 0, 
and models with (ii) Wq > —1 and Wi < have non-monotonic dark energy density px{z)- 
Models of type (i) have p'{z) < for 2 < Zcru and p'{z) > for 2 > Zcrit where 

Zcrit^ \{'i-+Wo)/wi\ . (10) 

Models of type (ii) have p'{z) > for 2 < Zcrit and p'{z) < for 2 > Zcru- 

All models other than types (i) and (ii) above have monotonic dark energy densities, 
decreasing with redshift for wq < —1 and Wi < and increasing with redshift for wq > —1 
and Wi > 0. 



IV. SIMULATED DATA 

We now construct simulated SN la data for dark energy Models 1, 2, 3, and 4 defined 
above, and investigate if we can recover the original theory from the simulated data. 

We simulate the data by distributing SNe la in z randomly per 0.1 redshift interval, with 
the total number per redshift interval as expected for SNAP. Here we assume that SNAP 
will obtain all SNe la in its survey fields up to 2 = 1.7 [25], similar to a supernova pencil 
beam survey [22,12]. We increase the number of SNe la at low redshifts, such that there 
are a minimum of 50 SNe la per 0.1 redshift interval at z < 0.5. We assume that these 
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additional low rcdshift sTipcrnovae will come from surveys of nearby SNe la. Thus each 
simulated data set consists of 2300 SNe la. 

The measured distance modulus for the l-th SN la is 

= + (11) 

where //^'^ = fip{zi) is the theoretical prediction in our dark energy model for a SN la at 
redshift zi [see Eq.(l)], and e*^'^ is the uncertainty in the measurement, including observational 
errors and intrinsic scatters in the SN la absolute magnitudes. In the simulated data set, 
we draw the dispersion e'-'-' for the l-th. SN la from a Gaussian distribution with variance 
^i^int = 0.16 mag. We simulate one set of data for each of the four models described in the 
Table. 

V. ESTIMATION OF DARK ENERGY FUNCTIONS 

We recover the dark energy function from each simulated data set. A number of tech- 
niques may be used to achieve this reconstruction. In this paper we use the adaptive iteration 
method ^ introduced by Wang & Garnavich (2001) [11] and Wang &: Lovelace (2001) [12] 
and described briefly here. Our current study builds on our previous work [5]. There we 
proposed a technique for determining the correct sign of p'x{z) if px{z) is not too close to a 
constant. To quantify how close to a constant we can go, in this paper we perform Monte 
Garlo simulations in order to obtain error bars for our results. 

A. Adaptive Iteration Method 

We start from the simulated data sets constructed from each of the four models described 
above. Given our data set, we now proceed as though we did not know which model it came 
from. We pretend we know nothing about the form of px{z). In attempting to reconstruct 
the dark energy density, we run through a series of test functions: we allow the test function 
Px^*(^) to be an arbitrary function. To approximate the function, we parametrize it by its 
value at + 1 equally spaced redshift values Zi {i = 0,1,2, . . . , N, zq = 0, zn — Zj^ax)- The 
values of p^x^{z) at other redshifts are given by linear interpolation, i.e., 

\Zi — Zi^l J \Zi — Zi^i J 

Zq — 0, Zn — Zynax- (12) 

The values of the dark energy density [i = 1,2,...,N) are the independent variables to 
be estimated from data. Again, we proceed as though we had absolutely no information on 
the function px{z), and treat it as a completely arbitrary function. 



^For a complementary method, see [24] which uses a Markov Chain Monte Carlo (MCMC) 
technique. 
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It is convenient to trade the + 1 parameters pi with the N parameters fi and the single 
parameter Q^, where 

fi^^ (i = l,2,...,iV) and = (13) 

PO Pcrit 



We define 



Po = Px(0) (14) 



and take Pcrit = SHq/ISttG) as the usual critical density. We thus have a total of N + 2 
parameters: the Hubble constant h = i/o/(100 km/s/Mpc), the matter energy density 
parameter Qm, and the N parameters /j (i = 1,2, . . . , A^) describing the test dark energy 
function. The complete set of parameters, then, is 

s = (h,nrn, Pi), (15) 

where i = 1, as described above. We will vary the number of bins between 1 and 14, 
and look for the optimal fit to the data. To illustrate, an arbitrary function may become a 
good approximation to the data for 4 bins whereas it is a miserable fit for 3 bins. 

In [5] we expanded the adaptive iteration method developed in Wang & Garnavich (2001) 
[11] and Wang & Lovelace (2001) [12] to include arbitrary time dependence of the dark energy 
density; unlike the earlier papers, we do not restrict ourselves to cases where p'x{z) > 0. 

The adaptive iteration method is designed to optimize the estimation of the dark energy 
density px{z) from data. It starts with fi = 1 for all i = 1,2,...,A^ (a cosmological 
constant), and builds fx{z) up iteratively while minimizing a modified statistics defined 
shortly. This adaptive iteration technique is further explained in the Appendix. 

We can now determine a best fit to the set of parameters s by using a statistic, with 

[8] 



X (s) = L -2 (16) 



where iJi^\zi\s) is the prediction for the distance modulus at redshift zi^ given the set of 
parameters s, and the sum extends over all the observed SNe la. Here ui is the dispersion 
of the measured distance modulus due to intrinsic and observational uncertainties in SN la 
peak luminosity. 

Assuming Gaussian errors, the probability density function for the parameters s is 



p(s) oc exp — — . (17) 




The normalized probability density function is obtained by dividing the above expression 
by its sum over all possible values of the parameters s. 

The probability density function of a given parameter Sj is obtained by integrating over 
all possible values of the other A^ + 1 parameters. To reduce the computation time, we can 
integrate over the Hubble constant h analytically, and define a modified statistic [11] as 
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X'^xl-^^{c, + hnl0^-2\nh\ (18) 
Here h* is a fiducial value of the dimensionless Hubble constant h, 

X^-E4W^^-/^?T' (19) 
I '^i 

^i-E4W^^-/^?0' (20) 

^2 = E;i> (21) 

with 

= ^Xp{zl■ h^h*)^ 42.384 - Slog /i* + Slog [HodL{zi)/c] , (22) 

The probability distribution function of the estimated parameters (excluding h) is now 
exp(— x^/2). It is straightforward to check that the derivative of with respect to h* is 

zero; hence our results are independent of the choice of h*. We take h* = 0.65. 

For a given choice of N, we can minimize the modified statistic of Eq.(18) to find the 
best fit and px{z) (parametrized by p^, i = 1, 2, N). 



B. Using the Value of Q,m to Constrain p'x{z) 

For the remainder of this Section, we restrict ourselves to Models 1, 2, and 3. We return 

to Model 4 in the Results Section below. 

To reiterate, we have started from three of the models defined in Section 111: (1) the cos- 
mological constant model with no time dependence in the energy density, (2) a quintessence 
model with p'x{z) > 0, and (3) an MP Cardassian model with p'x{z) < 0. We have con- 
structed a simulated data set for each of these models, and aim to see how well we can 
go backwards to determine the sign of p'{z) from this fake data set. In other words, can 
we reconstruct correctly the sign of the timc-dcpcndcncc of the energy density of the true 
model? To do this, we find (via the adaptive iteration technique) the model that best fits 
the observed matter density flm- 

For a given number of redshift bins N, we can minimize the modified statistic of 
Eq. (18) to find the best fit values fl^ and /f*^* for Qm and fx {z) = px (z) / px (0) parametrized 
by fi {i = 1, 2, N). For each model in Table 1, we obtain three sets of best fit parameters. 
We apply three different constraints to the arbitrary function p^x^^{z) in order to discover 
which one allows a good fit. The three constraints are: 
(i) Px^^{z) — Px^*(0) = constant; i.e., a cosmological constant model; 
(n) dp*^'\z)/dz > 0; 
and (ni) dp^^'\z)/dz < 0. 

We note that the second and third test functions are monotonically increasing and decreasing 
respectively, but do allow portions of the function to be flat as a function of redshift (i.e., 
constant for some but not all redshifts). 
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For each of these three constraints, we find the best fit parameters. If our technique 
works, the trial that gives the closest to the true corresponds to the correct un- 
derlying theoretical model (Model 1,2, or 3). For example, for the case where Model 2 
(p'xi^) ^ 0) is the theoretical model, if the trial with dp^x^{z) / dz > obtains the best value 
of then we have reproduced the correct time dependence of the dark energy density. 
Indeed we find that the technique works. 

Figure 1 shows our results: panels (a), (b), and (c) correspond to Models 1, 2, and 3. For 
each of the three models, the figure shows the best fit ^2^'*, under all of the three constraints 
above, for N values ranging from 1 to 14. The different constraints are represented by 
different point types. The dot-dashed horizontal line is our fiducial value of flm — 0.3 
(i.e., we are assuming that that this is the true value of the matter density), and the solid 
horizontal lines indicate 10% error bars about his fiducial value. We are assuming that flm 
is known to within 10% from other data sets. 

These plots are NOT intended to emphasize the dependence of OJ^^ on N, the total 
number of redshifts sampled via linear interpolation (see Eq. (12)). Indeed, as discussed 
above, the reason that we have found the best fit foi' a variety of N values is simply that 
the parametrization of the arbitrary function px{z) may be poor for one value of but 
excellent for another. For example, two points (A^ = 1) are perfectly adequate to describe a 
straight line function, but more points are needed to describe any more complicated function. 
Any one value of N may (by bad luck) give a bad result. We take a given model to be a 
good one if it lies within the 10% range on Vt^n for several values of N. The optimal value 
of is the one with the lowest x^, but we find that does not change much over a wide 
range of possible values of A^. We look for stability, i.e. for values of that do not change 
much as we vary A^ slightly. 

For Model 2 (with p'x{z) > 0) and Model 3 (with p'x{z) < 0), as A^ increases, the 
estimated values assuming the correct sign of p'xi^) asymptote to the "true" value of 
Vim = 0.3, while the estimated values il^f assuming the wrong sign of p'xi^) asymptote to 
an incorrect value of VLrn easily ruled out by observational constraints on VL^. Indeed our 
technique works. 

For the cosmological constant model (Model 1, top panel of Fig.l), the estimated Jl^* 
values assuming opposite signs for p'x{z) fall roughly symmetrically on opposite sides of the 
"true" value, and they lie within the error bars on Q„i. This indicates a degeneracy between 
A models and dark energy models that have mildly time-dependent dark energy density. In 
other words, it will be difficult to differentiate a constant dark energy density from one that 
has a very small dependence on redshift. We will examine this possible degeneracy further 
in Section VI to see how well one can differentiate between these two alternatives. 

In conclusion, we can determine the correct sign of p'x{z) if Px{z) is not too close to a 
constant. To quantify how close to a constant we can go, we need to add error bars to the 
points in figure 1. 

C. Using Monte Carlo to Determine Errors 

We evaluate errors by simulating random fluctuations around a fiducial model. A possible 
choice of fiducial model would be the input (theoretical) model we adopted to generate the 
SN la data set in the first place. However, to be closer to a realistic situation in which 
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the underlying model has to be determined from data, we will choose our fiducial model 
to be the model that best fits the data (for us, the simulated data). Specifically, this 
is the cosmological model with p^{z) (BF = best fit) given by the best-fit values pf^ 
{i — 0, 1, ...,iV) determined above. We have described in the previous section the method 
by which we obtain our best fit dark energy densities for each of Models 1, 2, and 3. Here, 
we add random errors to the "measured" distance moduli. 

To derive robust error distributions of the estimated parameters Q'^^^ and f^^^ = Pi/po 
{i = 1, 2, A^, see Eq.(12)) from each data set, we create 10"^ Monte Carlo samples by 
adding dispersion in peak luminosity of Arriint = 0.16 mag to the distance modulus Pp{z) 
[see Eq.(l)] predicted by the best-fit model (i.e., assuming that the best-fit model is the true 
model). This is equivalent to making 10^ new "observations", each similar to the original 
data set [21]. The same analysis used to obtain the best-fit model from the data is performed 
on each Monte Carlo sample. The distributions of the resultant estimates of the parameters 
(O^* and f^^^) can be used to derive the mean and confidence level intervals of the estimated 
parameters. Wang & Lovelace [12] showed that such a Monte Carlo analysis gives less biased 
estimates of parameters than a maximum likelihood analysis, i.e., the Monte Carlo mean of 
estimated parameters deviate less from the true values of the parameters. 

VI. RESULTS 
A. Models 1, 2, and 3 

Starting from the three dark energy models 1, 2, and 3 in Table 1 with p'xi^) — 0; 
p'x{^) — ^'^d p'x{z) < respectively, we obtained simulated data sets. Using these 
simulated data sets^, we worked backwards to estimate the best fit values of the + 1 
parameters as defined in the previous section. In particular, for each of the simulated data 
sets we made three trial assumptions, i.e., constant, increasing, and decreasing dark energy 
density. For each assumption we found an estimate of Vim and then ran 1000 Monte Carlo 
samples to obtain error bars. The results of this analysis are shown in Figure 2. The allowed 
value of VLmi assumed to be known to within 10% is shown with arrows: VL^n = 0.3 ± 0.03. 

Fig. 2 confirms the conclusions that we made based on Fig.l. If the true dark energy 
density varies with time (quickly enough) and monotonically (see Figs. 2b and 2c), the dark 
energy models with the wrong sign of p'x{z) or with constant px{z) are easily ruled out at 
> 95% C.L. if f2„j is known to ~ 10% accuracy. If the true dark energy density is constant 
with time (Fig. 2a), then while the correct model (A model) is favored, the incorrect models 
(with px{z) either increasing or decreasing with redshift) could imitate the correct model 
if the time variation in px {z) is small enough. We will quantify the size of the degeneracy 
region shortly (see Fig. 5 and related discussion). 



^Note that here, once data become available, the "simulated data set" will be replaced by the real 
data set. We consider three different simulated data sets based on different dark energy models, 
since the nature of dark energy is unknown. 
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Note that due to computational constraints, we have chosen N = 6 for demonstration. 
For an actual data set, the Monte Carlo analysis outlined here should be performed for all 
values of N (from 1 to a reasonably large number) . 

Figure 3 shows the best fit dark energy density for the simulated data set based on 
Model 2 (which has p'x{z) > 0). To obtain figure 3, we take advantage of the fact that our 
analysis in Fig. 2(b) has allowed us to extract the sign of p'xiz)-, with the assumption that we 
know to 10% accuracy. Once the sign of the time dependence is known, we can deduce 
the best fit dark energy density px{z)/ px{0) shown in Fig. 3, estimated from the Monte 
Carlo analysis of the simulated data set based on Model 2. The solid hue is the true model, 
i.e., the theoretical curve. One can see that the dotted points with error bars, obtained 
from the simulated data using our technique, match the true model very well. The value 
of the matter density corresponding to Fig.3 is Qm = 0.314(0.271,0.341) (mean and 68.3% 
confidence range). The corresponding figure for Model 3 {p'x{z) < 0) was published in [5]. 
Again, one can see that the time dependence of the energy density can be determined quite 
well. 

B. Model 4 (linear equation of state) 

Lastly we come to theoretical dark energy models with a linear equation of state Wx{z) = 
Wo + Wiz. Again, we simulate data based on this underlying model, and again we ask how 
well we can determine the time dependence of the energy density of each model. This linear 
form of the equation of state is the most common parametrization used by the community. 
As stressed previously, one can extract far more information from data if we aim to learn 
about the properties of the underlying energy density, in this case its time dependence. In 
fact, because of our ignorance of the true nature of dark energy, it is dangerous to rely on 
specific parametrizations. Hence, we perform the analysis described in the paper, in which 
we treat the dark energy density px{z) as a free function in extracting information from the 
data (in this case simulated data) . We will do a blind test (non-parametric study) , in which 
we do not assume the linear form of the equation of state, to see whether the time-variation 
of dark energy density can be ascertained. This will establish a point of reference with the 
work by others. In other words, we generate models with a given equation of state but then 
reconstruct (the time dependence of) the dark energy density itself. 

We studied a grid of models, with —1.2 < Wq < —0.5, and —1.5 < Wi < 0.5, and a 
grid spacing of = 0.1 (0.2 for Wq < —1) and Awi = 0.1. This represents a total of 
147 models. For each of these theoretical models, we followed the same procedure described 
above. First, we simulated data based on each model; next we obtained the best fit p'x^{z) 
for each model with the adaptive iteration method (the best fit function chosen due to 
its producing Q^* closest to the correct value); finally we created 10^ Monte Carlo samples 
(fluctuations about the the best fit function) to obtain error bars. We followed this procedure 
for each of the 147 models. Fig. 4 shows the estimated fi^f values (with la standard 
deviations) for a subset of the wx{z) = wq + wiz models that we have studied, assuming: 
(1) p'^'^iz) = constant; (2) dp'^'^{z)/dz > and (3) dp^^^^{z)/dz < 0. The dotted lines in 
each figure denote flm — (0.291,0.309), i.e., assuming that flrn is known to 3% from other 
observations. 

Fig. 4 runs through a subset of the grid of theoretical models that we have studied. 
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namely the models with (a) wq = —1.2, (b) wq = —1, and (c) wq = —0.8. In each plot we 
have selected one value of Wq and show results for a variety of values of Wi. Again, these 
values of Wq and Wi correspond to the underlying theoretical model. Based on each of these 
sets of Wo,Wi we simulate data and show how well we can ascertain the time dependence 
of the dark energy density. We plot the value of with error bars that results from 
our different trial assumptions, as a function of different (theoretical) values of wi in the 
underlying models. We have obtained three different i^^^'s (for the three different trial 
assumptions) for each underlying Wi. 

In Fig. 4a, we have taken the underlying set of theoretical models to have wx{z) = 
—1.2 + wiZ. In other words we took wq — —1.2 and allowed Wi to vary from —1.5 to 0.5 
with spacing Awi = 0.1. Our goal is to obtain the correct value of only for the test 
function that has the correct sign of p'x{z)- We see that our technique has indeed reproduced 
the correct answer, p'x{z) < for wi < 0; however, the answer is ambiguous for wi > 0, 
where dark energies that are constant or increasing with redshift seem to give equally good 
fits. In fact we understand the reason for the success of the technique for some values of 
Wi and the failure in other regimes, and explain it here. As a reminder, we showed below 
Eq.(9) that p'x{z) has the same sign as 1 + Wq + Wi{z). Hence, for wq = —1.2, the sign 
of p'x{z) is equal to —0.2 + Wiz. As we discussed in the last paragraph of Section III, for 
Wo < (which is the case for Fig. 4a), the dark energy density is monotonic with redshift 
for wi < 0, but nonmonotonic for Wi > 0. However, as our trial functions we have only 
considered monotonic dark energy densities. Hence it is not surprising that our technique 
only reproduces the correct time dependence of the dark energy in the regimes of parameter 
space that correspond to dark energies that are monotonically increasing or decreasing in 
time. The technique works within the regime of validity of our trial assumptions. 

In Fig. 4b, the underlying model is wxiz) = — 1 + WiZ, again for wi ranging from —1.5 
to 0.5. In this case the corresponding dark energy density is monotonic with redshift for all 
values of Wi, and the technique reproduces the right answer for all Wi. 

In Fig. 4c, the underlying model is wx{z) = —0.8 + WiZ, again for Wi ranging from 
— 1.5 to 0.5. This corresponds to a dark energy density that is monotonic for Wi > but 
non- monotonic for Wi < 0. The dark energy density grows with redshift until Zcru = (1 + 
Wq)/\wi\ = 0.2/|wi| (as obtained from Eq.(lO)) and then decreases with redshift. However, 
the value of Zcj-u is so small in most models as to be irrelevant; i.e. the function is essentially 
monotonically decreasing over most of the redshifts at which supernova data are taken. For 
example, for wi = —1.5, Zcru = 0.133. As Wi approaches 0, the value of Zcru grows; clearly 
at Wi — —0.2 and Zcrit — 1, the turnover from a growing dark energy density to a decreasing 
one takes place in the middle redshift of the supernova data. In that case we would expect 
a constant trial with Pxi^) = to be no worse a fit than a monotonically increasing or 
decreasing one. In other words, a cosmological constant appears to fit the data, although of 
course it is the wrong answer. Indeed this is borne out by the results of Fig. 4c. Here, our 
technique is very successful at reproducing the correct sign of the time dependence of the 
dark energy in those regimes where it is sensible to approximate the dark energy density as 
being monotonic in time. 

Models with Wo = —0.9, —0.7, —0.6, and — 0.5 behave similarly to those shown in Fig.4. 
As Wo gets less negative, the value of Zcrit (at which the derivative of the dark energy density 
changes from positive to negative) gets bigger. Again, our technique works to obtain the 
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correct time dependence of the dark energy density if it is monotonic in time. 

To quantify how much the p'x{z) 7^ models differ from the A models, we have defined a 
quantity Q as the number of standard deviations in the difference of the average estimated 
value of ftrn with constant and non-constant dark energy density, 



Q 







pji:=non— const ^ 


\ PX-- 


=const ^ 






Px =non— const ) 


' + (AO-* 


y 

px=const/ 



(23) 



The thick solid line in the bottom panel of each section of Fig. 4 shows Q as a function of 
wi for various wq values. As there are two test functions with time varying dark energy 
density, p'x{z) > and p'x{z) < 0, in the figures we plot the quantity Q for that time 
varying function that produces the value of that best fits the data. 

We note that our work is quite different to prior studies performed by other authors. 
Previous work (e.g., [9,10,26,27]) has attempted to examine how well one can reconstruct the 
dark energy equation of state, particularly if one assumes it has the form wxiz) — Wq + WiZ. 
Our study differs from previous studies in being non-parametric (hence not dependent on 
Wxiz) = Wq + Wiz being the correct model, but only using it as a test case), and spanning a 
continuous parameter space in (wq, wi) (via smooth interpolations in our grid of 147 models). 



How Small Can We Make the Region of Degeneracy? 

We can ask the question: what is the range of {wq, wi) over which one cannot determine 
the sign of the time dependence? In other words, what is the range of degeneracy between a 
cosmological constant and a time-changing vacuum energy appearing to fit the data equally 
well? And how far can it be shrunk down? 

Fig. 5 shows the [wq, Wi) parameter space that we have studied. The models that lie 
within the shaded region cannot be differentiated from a A model even if O^ is known to 1%. 
Similarly, models which lie within the dotted, dashed, and solid lines cannot be differentiated 
from a cosmological constant model if O^ is known to within 3%, 5%, and 10% accuracies 
respectively. 

The degeneracy region shown in Fig. 5 is centered about the line 

l + wo^ 3~^i' 

where z^ax is the maximum redshift of the survey {z^ax = 1-7 for SNAP). 

In the parameter space outside of the shaded region in Fig. 5, p'x{z) 7^ models are 
preferred over the A models if O^ is known to 1% accuracy, indicating our ability to detect 
the time-variation of the dark energy density at lu or higher significance levels. 

Finally, we point out a way to reduce the degeneracy by taking advantage of the fact 
that the slope in Eq.(24) changes for different values of Zmax- Hence changing z„iax will 
rotate the degeneracy region in {wq^Wi). One can choose a variety of different vahies of 
Zmax to break the degeneracy. In other words, in addition to using the entire data set, one 
can restrict the data out to a variety of different cutoff redshifts to obtain complementary 
information. If, in addition to the full data set, we consider only those data to a second 
maximum redshift Z2,max-i we can reduce the degeneracy region in Fig. 5 significantly. In 
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Figure 5 we have plotted both the degeneracy region obtained using all the data, as well as 
dot-dashed lines illustrating the different degeneracy region if only data out to Z2,max = 0.5 
were used.'^ The combination of information from these dot-dashed lines together with the 
shaded region allow us to break the degeneracy substantially. The shaded region in Fig. 5 
bounded by dot-dashed lines illustrates the smallest possible degeneracy region if is 
known to 1%. 

VII. CONCLUSIONS 

We have investigated how well future supernova data from dedicated telescopes (such 
as SNAP), when combined with data of nearby supernovae, can be used to determine the 
time dependence of the dark energy density. For quantitative comparison, we have done 
an extensive study of a number of dark energy models, with dark energy density that is 
constant, increasing, and decreasing in time. Based on these models we have simulated data 
sets in order to show that we can indeed reconstruct the correct sign of the time dependence 
of the dark energy density. 

Among the dark energy models we studied are those parametrized by an equation of state 
wx{z) = wq + wiz. Here, w = p/ p. We studied a grid of 147 models, for —1.2 <Wo< —0.5, 
and —1.5 < Wi < 0.5. We emphasize that it is the dark energy density that we reconstructed, 
not the equation of state. We find that there is a degeneracy region in the {wo, Wi) parameter 
space centered near 1 + wq = —wiZmax/'^ (where Zmax is the maximum redshift of the survey, 
e.g., Zmax — 1-7 for SNAP); the models that lie within this region cannot be differentiated 
from a A model even if 0,^. is known independently to 1% accuracy (we compute the size of 
the region for VL„i known to varying degrees of accuracy). Outside of this degeneracy region, 
we can detect the time variation of the dark energy density at Icr or higher significance 
levels. 

We emphasize that, given the same data, one can learn much more by reconstructing the 
dark energy density directly (and its time dependence) than by attempting to reconstruct 
its equation of state. 
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APPENDIX: ADAPTIVE ITERATION TECHNIQUE 



^These are qualitative illustrations, not from actual calculations (which would involve lengthy 
computations) . 
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The goal of the adaptive iteration technique is to reconstruct (as accurately as possible) 
the function fxiz) = px{z)/ px{0) from a simulated data set. We start from a time inde- 
pendent function (in which fx{z) = 1 for all z) and build up the function iteratively to find 
that function which best matches the data. 

To illustrate the adaptive iteration technique, we present an example. In this appendix 
we here restrict our discussion to monotonically increasing forms of fx{z). Let us consider 
the case where we break up the function (and the simulated data) into five equally spaced 
bins in redshift space, i.e., = 5 so that z ranges from 0, Zi, Z2, z^, Z4, z^. We will start with 
a flat function, fx{z) = px{z)/ px{0) — 1 and build up from there. In all the iterations, 
we will always keep /x(0) = Pxi^ — 0)/px{0) — 1 fixed (it's an identity equation since 
px{0) = px{z — 0)), and vary fx{z) at the other values of z. Here is how we proceed in the 
first iteration: 

(1) Compute for fx{z) = 1- 

(2) Increase fx{zi) for all i — 1,2,3,4,5 by one stepsize. A, i.e., fx{zi) = 1 + A for all 
i — 1,2,3,4,5. Compute x^. If the current is smaller than the previous (from the 
previous step), the new fx{zi) values are favored; keep them. Otherwise, the previous values 
of fx{zi) = 1 are favored. As stepsize we used primarily A = 0.01, but also A = 0.05 and 
0.1 for comparison. We find that the results are not sensitive to the stepsize; of course, the 
smaller the stepsize, the longer the running time. 

(3) Increase fx{zi), i = 2, 3, 4, 5 by one stepsize. A, i.e., fx{zi) = fx{zi)P^''''+A, i = 2, 3, 4, 5, 
where fxizi)^^^^ are the fx{zi) values favored by the previous step. Compute x^- If the 
current is smaller than the previous Xmin (from the previous step), the new fx{zi) values 
are favored; keep them. Otherwise, fxiziY^^^ are favored. 

(4) Increase fx{zi), i = 3,4, 5 by one stepsize. A, i.e., fx{zi) = fx{ziY''''" + A, i = 3,4,5, 
where fx{zi)^^'" are the fx{zi) values favored by the previous step. Compute x^- If the 
current x^ is smaller than the previous Xmin (from the previous step), the new fx{zi) values 
are favored; keep them. Otherwise, fxiziY"^^"" are favored. 

(5) Increase fx{zi), ? = 4, 5 by one stepsize, A, i.e., fx{zi) = fxiziY^'^^ + A, -i = 4, 5, where 
fx^ZiY'^^'" are the fx{zi) values favored by the previous step. Compute x^- If the current 
X^ is smaller than the previous Xmin (fi'om the previous step), the new fx{zi) values are 
favored; keep them. Otherwise, fxiz^y^'^'" are favored. 

(6) Increase fx{zi), i = 5 by one stepsize. A, i.e., fx{z.i) = fxiziY^'^^ + A, i = 5, where 
fx{ziY^^^ are the fx{zi) values favored by the previous step. Compute x^- If the current 
X^ is smaller than the previous Xmin (from the previous step), the new fx{zi) values are 
favored; keep them. Otherwise, fx{zi)'^^'" are favored. 

In the second iteration, repeat steps (2)- (6), but with these changes in step 2: replace 
fx{zi) = 1 + A with fx{zi) = fxiziY"''^'" + A, replace "the previous x^" with "the previous 
Xmin" J and replace "the previous values of fx{zi) = 1" with fx{ziY'''^^'\ 

Subsequent interations follow the same procedure as the second iteration. Continue the 
interations until Xmm stops changing. 

We successively perform further iterations to ascertain the function px{z) that best fits 
the data. As described here, one can only end up with a monotonically increasing form 
for p{z). For monotically decreasing function, the procedure is exactly the same, with -|-A 
replaced by —A. 
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FIG. 1. The estimated best fit shown for a variety of values of the number of parameters 
N. The three panels arc for the dark energy models 1, 2, and 3 (given in Table 1) respectively 
Each panel plots the estimated il^^* for three trial functions: crosses indicate constant /9^*'*(z), 
squares indicate dp^^^ /dz > 0, and triangles indicate dp^^^ /dz < 0. The dot-dashed line indicates 
the "true" value of = 0-3 in each model, the solid lines indicate the ±10% range of In 
each panel the test function with the correct time dependence (same as the underlying theoretical 
model) produces an acceptable fi^* that matches data; the test function with the wrong time 
dependence produces an incorrect value of J7„i and hence can be ruled out. This technique can 
reproduce the correct time dependence of the dark energy density. 
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FIG. 2. (a) Distributions of the estimated fi^ values from 10^ Monte Carlo samples. Here, the 
underlying theoretical model is a cosmological constant. The Monte Carlo samples are obtained 
under three different trial assumptions about the time dependence of the dark energy p'xi^)'- 
increasing, decreasing, or constant. The allowed value of = 0.3 it 0.03 is indicated. The correct 
value of is reproduced for the correct time dependence of the dark energy, p'j^ = 0; hence one 
is led to conclude that the underlying model is a cosmological constant. Note, however, that for 
this case some of the models with time dependent px also produce the right values of Qm, so that 
there is some degeneracy. 



19 



500 



T r 



(b) 



allowed Q. 



T 1 1 1 1 1 1 r 

Theoretical Model 2: 



400 hsNAP 

10^ Monte Carlo samples 

Px'(z)=0 

Px'(z)^0 
Px*(z)^0 



^ 300 



n 200 



correct model 



100 - 







-r 1 

I t 1 

I I L 

1 1 1 r 

I.I I 1 !-in 



0.2 



0.3 



0.4 



0.5 



m 



FIG. 2. b) Same as Fig. 2(a), but for the simulated data set based on Model 2. Upon running 
10^ Monte Carlo samples, the correct value of Q„i, is reproduced only for the trial functions with 
p'x > 0, which agrees with the time dependence of the underlying model. [The underlying model 
is a quintessence model with Wx{z) = — 1 + 0.5z and p'x{z) > 0.] We have taken A?^ = 6 for 
illustration. Hence one recovers the correct time dependence of the energy density. 
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FIG. 2. c) Same as Fig. 2(a), but for the simulated data set based on Model 3. Upon running 
10^ Monte Carlo samples, the correct value of Q„,, is reproduced only for the trial functions with 
p'x < 0, which agrees with the time dependence of the underlying model. [The underlying model 
is a MP Cardassian model with n=0.2 and q=2, and p'x{z) < 0.] We have taken N = 6. Hence 
one again recovers the correct time dependence of the energy density. 
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FIG. 3. Best fit dark energy density px{z) / px{^) estimated from the Monte Carlo analysis of 
the simulated data set based on Model 2, which has p'xi^) > 0- This plot assumes that il„,, is 
known to 10% accuracy [and that the sign of the time dependence has been extracted as shown 
in Figure 2]. The soHd line is the true model, i.e., the theoretical curve. One can see that the 
dotted points with error bars, obtained from the simulated data using our technique, match the 
true model very well. Here, we have taken N = Q. 
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FIG. 4. The estimated Qm values (with la standard deviations) for wx{z) = wo + wiz models. 
In each plot in panels a-g, we have selected one value of wq and show results for a variety of 
values of wi; these values of wq and wi correspond to the underlying theoretical model. Based 
on each of these sets of wo,wi we simulated data, and in the upper half of each panel, plotted 
the values of with error bars that result from our three different trial assumptions, as a 
function of different (theoretical) values of wi in the Tindcrlying models. The three different trial 
assumptions are represented by: circles for p^^* = constant, squares for dp^^^ /dz > 0, and triangles 
for dp'^^/dz < 0. The dotted lines denote 0„ = (0.291, 0.309) (here, we assume that O^, is known 
to 3% from other observations.) In the lower half of each panel, we have plotted the quantity Q 
defined in Eq.(23), the number of standard deviations in the difference of the average estimated 
vahic of with constant and non-constant dark energy density. 

(a) Underlying theory: wq = —1.2 and —1.5 < Wi < 0.5. The correct value of $7^ is obtained for 
p'xi^) — ^1 < 0; here our technique has indeed reproduced the correct time dependence of 
the dark energy. However, the answer is ambiguous for wi > 0, where the underlying theory is 
non- monotonic . 
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FIG. 4. (b) Underlying theory: ivq = —1 and —1.5 <wi< 0.5. Here our technique obtains an 
allowed value of and hence reproduces the correct time dependence of the dark energy for all 

Wl. 
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FIG. 4. (c) Underlying theoretical model: wq = —0.8 and —1.5 < u^i < 0.5 (see discussion in 
text) . 
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FIG. 5. The (wq, wi) parameter space that we have studied. For the theoretical models with 
wq and wi that lie within the shaded region, the reconstructed p'xi^) cannot be differentiated from 
a A model (indicated by a fat circle) even if Q^n is known to 1%. Outside of this region, the sign 
of p'x{z) can clearly be ascertained (if rim is known to 1%). Similarly, models which lie within the 
dotted, dashed, and solid lines cannot be differentiated from a constant p'x{z) = model if Q,m 
is known to within 3%, 5%, and 10% accuracies respectively. The degeneracy region is centered 
about the line 1 + it;o — — ^"^""^ wi, where Zmax is the maximum redshift of the survey {zmax = 1-7 
for SNAP). The dot-dashed line illustrates (roughly) the different degeneracy region if only those 
data out to a cutoff redshift of 0.5 were used, if were known to 1%. Note that the degeneracy 
can be reduced by examining different portions of the data out to different redshifts. 
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